Weak Local Adaptation to Climate in Seedlings of a Deciduous Conifer Suggests Limited Benefits and Risks of Assisted Gene Flow

ABSTRACT Assisted migration provides a potential solution to mitigate the increasing risks of forest maladaptation under climate change. Western larch (Larix occidentalis Nutt.) is a deciduous conifer species undergoing assisted migration beyond its natural range in British Columbia into areas that have become suitable based on climatic niche modelling. We established a seedling common garden experiment in raised beds in a warm location outside the natural range for three growing seasons, with 52 natural populations from across the species range and 28 selectively bred families from British Columbia. Intraspecific genetic variation in growth, phenology and cold hardiness was analyzed to test for signals of local adaptation and the effects of selective breeding to better understand the implications for assisted migration and breeding for future climates. We found weak differentiation among populations in all traits, with the proportion of additive genetic variance (Q ST) ranging from 0.10 to 0.28. Cold hardiness had the weakest population differentiation and exhibited no clines with geographic or climatic variables. Selective breeding for faster growth has maintained genetic variation in bud flush phenology and cold hardiness despite delaying bud set. The weak signals of local adaptation we found in western larch seedlings highlights that assisted gene flow among populations is likely to have limited benefits and risks for mitigating maladaptation with climate change. Our findings suggest that assisted migration outside of the range and selective breeding may be important management strategies for western larch for future climates.


| Introduction
Forest trees are becoming increasingly vulnerable to climatic stress with rapid climate change (Allen et al. 2015;Hartmann et al. 2022).Shifting climates and greater climatic variation are disrupting historical patterns of local adaptation in tree populations, increasing the risk of maladaptation and the cascading ecosystem-level consequences of declining forest health, and presenting new and urgent challenges for forest management (Millar and Stephenson 2015).A potential management solution for mitigating the increasing risk of maladaptation is assisted migration, the managed translocation of genetic material to areas where it is forecasted to be optimally adapted in the future to facilitate more rapid adaptation to climate change (Aitken and Whitlock 2013).Assisted migration can include the movement of genetic material among populations within the current range, that is, assisted gene flow, and assisted migration beyond the current range (Aitken and Whitlock 2013).Tree populations have limited capacity to adapt or migrate quickly enough to track current and future climate change because of their long generation times but they also tend to have high gene flow and limited population isolation, making them well suited for assisted migration.
To evaluate the potential impacts of assisted migration, it is critical to understand patterns of local adaptation to historic climates in trees.More than 250 years of common garden experiments in forestry have shown that local adaptation to climate is common in tree species and that the strength and scale of local adaptation varies among species (Alberto et al. 2013; Leites and Benito Garzón 2023).In common garden studies, population differentiation in relevant traits and the association between traits and source environments is regarded as evidence consistent with local adaptation to climate (Aitken and Bemmels 2016;Wadgymar et al. 2022;Leites and Benito Garzón 2023).Many widely distributed and economically important conifer species exhibit strong signals of local adaptation to climate, while many broadleaf and some conifer species show weak to no signal, suggesting a relatively strong role of phenotypic plasticity to accommodate environmental variation (Leites and Benito Garzón 2023).Phenotypic plasticity is the capacity of a genotype to produce different phenotypes in response to environmental variation.Adaptive phenotypic plasticity can help stabilize fitness traits and the limits of phenotypic plasticity can determine how much environmental change a genotype can withstand before fitness is compromised (Ghalambor et al. 2007).Phenotypic plasticity can buffer populations against environmental change, allowing persistence in the short term, but it can delay genetic adaptation in the long term (Valladares et al. 2014).Understanding the strength and scale of local adaptation relative to phenotypic plasticity in different tree species can help inform effective strategies for assisted migration in future climates.
It is well established that cold temperatures are important drivers of local adaptation in temperate trees, especially conifers, which have a long evolutionary history of adaptation to harsh climatic conditions that has allowed them to dominate temperate and boreal forests throughout the world (Bond 1989).Climate change and assisted migration approaches are likely to shift the timing and distribution of late spring and early fall frost risks relative to the seasonal developmental cycle of locally adapted tree populations (Sang et al. 2021).One of the major challenges of assisted migration for temperate and boreal tree populations is the need to balance the potential risks of frost damage and compromised survival at early seedling stages against the benefit of mature trees optimally adapted to warmer climates decades to centuries later (O'Neill et al. 2014;Erlichman et al. 2024).Further, cold temperatures are important cues for the process of cold hardening in trees in the fall, and warming or more variable temperatures could delay this process resulting in greater vulnerability to cold injury (Bansal et al. 2015).Understanding the potential impacts of climate change and seed transfer on local adaptation to cold temperatures and seasonal frost risks is therefore critical for minimizing risks associated with assisted migration.
Locally adapted temperate and boreal tree populations have evolved timing of active growth so that trade-offs between growth and survival due to seasonal frost risks are optimized.While trade-offs between growth and seasonal frost tolerance are well established in natural populations (Howe et al. 2003), trade-offs or correlated responses associated with selective tree breeding for faster growth are less clear yet the risks associated with compromised frost tolerance are an increasing concern for future tree breeding and assisted migration under rapidly changing climates.Traditional tree breeding programs aim to increase growth gains by selecting for genotypes with faster growth rates, and this may result in an extended growing period.When selection delays growth cessation, there is a risk of a phenological mismatch with seasonal frost risks, resulting in selectively bred trees that are more vulnerable to cold injury compared to natural populations.Alternatively, if selective breeding enhances growth without strong impacts on phenology or cold hardiness, the synchrony with seasonal frost risks could be maintained and seed transfer based on climate could be similar between natural and selectively bred sources.Encouragingly, previous common garden studies of lodgepole pine and interior spruce have found that selective breeding for faster growth has mostly maintained variation in cold hardiness (MacLachlan et al. 2017;MacLachlan, Yeaman, and Aitken 2018).
We studied genetic variation in traits related to the seasonal development cycle of western larch (Larix occidentalis Nutt.), a deciduous conifer species undergoing assisted migration trials beyond its current natural range in British Columbia (Rehfeldt and Jaquish 2010;O'Neill et al. 2017).Climate niche modelling has projected that areas of suitable climate for western larch have expanded beyond and contracted within its current range and will continue this trend in the coming decades, making it a desirable candidate for assisted migration beyond the range (Rehfeldt and Jaquish 2010;MacKenzie and Mahony 2021).We established a seedling common garden experiment in raised beds in a location warmer than the natural range for three growing seasons with a range-wide collection of 52 natural populations to simulate effects of assisted migration and warming.We included 28 selectively bred full-sib families to compare with natural populations from the same regions to test for the effects of selective breeding on climate-relevant traits.We measured traits related to the seasonal developmental cycle (height growth, bud phenology, and cold hardiness) to address the following questions: (1) How much population differentiation exists in climate-relevant traits among natural populations?(2) Does trait differentiation among populations exhibit clinal variation with climatic or geographic gradients?(3) How does selective breeding for faster growth affect variation in climate-relevant traits?

| Study Design
We obtained 52 open-pollinated seedlots from natural populations across most of the natural range (Figure 1).We selected populations from available seedlots dispersed in geography and climate spanning a range of 6° latitude, 7° longitude, 1100 m elevation, 6°C mean coldest month temperature, and 1325 mm mean annual precipitation (Table S1).The common garden site was situated outside of the natural range in Vancouver, BC, which has a warmer climate (Figure S1).Seedlots were obtained from stored collections from the Rocky Mountain Research Station, Inland Empire Tree Improvement Cooperative, and British Columbia Tree Seed Centre.Each seedlot was bulked from at least 10 maternal trees in British Columbia and 5-10 maternal trees in the United States, generally sampled at least 100 m apart at a given site.Populations of the interior portion of the range (i.e., Rocky Mountain region) are well represented but populations from the westernmost portion of the range (i.e., Washington and Oregon) were not included in the study due to the small number of maternal trees sampled per location for the few available seedlots.To test for the effects of selective breeding on traits relative to natural populations, we obtained seed from 28 selected full-sib families of the western larch breeding program of the British Columbia Ministry of Forests from two breeding zones established in 1990 (Figure S2).Selected full-sib families included 32 parents total, 28 of which were represented in more than one family, allowing for estimation of additive genetic variation and narrow-sense heritability values.
Seeds were soaked in water and cold stratified for 3 weeks, then sown in 55 cm 3 Stuewe & Sons Ray Leach UV containers filled with Sunshine Mix #4 potting media on April 13-16, 2019.Two to four seeds were sown per container depending on expected germination rates of each seedlot and thinned to one seedling after germination.Seedlings were grown in a greenhouse under 18-h days and an average temperature of 25°C.Trays of 200 seedlings each were kept well-watered, fertilized with Peters Excel 11-41-8 N-P-K water soluble fertilizer, and rotated weekly to avoid any spatial effects inside the greenhouse.In order to induce bud set and dormancy, all seedlings were put under a 10-h photoperiod from September 30th onward.By November 13th, all seedlings had set bud and were moved to controlled chambers set at 4°C in the dark to simulate winter and chill seedlings to break bud dormancy.
We transplanted seedlings to raised outdoor beds in April 2020.Seedlings were planted in four blocks of 240 plants in a complete randomized block design.Seedlings were planted with 4 cm spacing, with three randomly assigned individuals per seed source within each block, 12 plants per seed source, totaling 960 seedlings.Each block was surrounded by a row of buffer seedlings that were not measured to minimize edge effects.
Seedlings were watered through the growing season and we applied Peters Excel 15-5-15 N-P-K water soluble fertilizer at the beginning of the growing season to avoid nutrient deficiency.Signs of a fungal pathogen infection, identified as Sirococcus conigenus, were detected in a small number of seedlings in early June 2021, causing the terminal shoot to curl downwards, the needles to turn brown and eventually fall off.A broad-spectrum fungicide treatment, Aliette, was administered immediately and plants were regularly monitored for signs of infection.No signs of infection spread were detected after treatment.Seedlings with dead tops resulting from the infection (<5%) were excluded from analyses.

| Climate and Geographic Data
We used the geographic coordinates of each population to obtain interpolated climate records for the period of 1961-1990 from ClimateNA (Wang et al. 2016).This is the earliest period that weather stations produced reliable data and was assumed to be early enough to represent historic conditions.For testing phenotype-environment associations, we selected five annual climate variables based on biological relevance and pruned for collinearity for phenotype-environment associations: mean coldest month temperature (MCMT), continentality (TD) calculated as the difference between mean warmest month temperature and mean coldest month temperature, frost-free period (FFP), mean annual precipitation (MAP), and summer-heat moisture index (SHM).We also analyzed trait associations with geographic variables of latitude, longitude, and elevation.

| Trait Data
We measured six traits related to the seasonal development cycle in conifers (height growth components, bud flush, bud set, and fall cold hardiness).Phenotypic data from the third year of growth ( 2021) was analyzed for this study.Height was measured five times throughout the growing season and final height was measured after growth cessation and final bud set on September 6, 2021.Early growth was estimated as the height on June 21, 2021 (the date seedlings began to set bud for the first time) minus the initial height measured at the beginning of the growing season on May 13, 2021.Early growth includes both determinate growth from shoot primordia that overwintered in buds and indeterminate free growth from neoformed leaves and internodes.In western larch, free growth extends beyond the juvenile stage and remains a significant component of growth potential in mature trees (Schmidt and Shearer 1990).We could not separate determinate and free growth in this study because preformed and neoformed leaves are not readily distinguishable based on morphological differences; therefore, we analyzed early and late season growth instead.Late growth was estimated as final height minus the height on June 21, 2021.Late growth includes free and lammas growth, a form of indeterminate shoot growth resulting from reflushing after initial bud set in the current growing season.
Visual assessments of the timing of terminal bud flush and bud set were collected to analyze the timing of initiation and cessation of height growth, respectively.A binary system was used to record the visual observations.Bud flush was determined as the parting of bud scales and the emergence of needles and bud set as the formation of visible brown bud scales on apical buds.Several seedlings set bud multiple times due to lammas growth, and final bud set was determined as the date after which no reflushing was observed.Observations of terminal bud phenology were recorded every 5 days at the beginning of the growing season until all seedlings had fully flushed from March 8 until April 13, 2021.Bud set was monitored weekly from June 21 to September 22, 2021.Bud flush and bud set dates were analyzed as day of the year starting January 1.
Cold hardiness was measured in October as minimum air temperatures began to approach 0°C at the test site.Measurements involved artificially freezing and measuring electrolytic leakage of branch tissue following a protocol adapted from (Hannerz et al. 1999).Three samples, each comprising three branch segments 5 mm length, were collected from each seedling; two samples were used for two different freeze test temperatures and the third sample was an unfrozen control.Test temperatures were selected based on preliminary freeze tests to result in approximately 50% cold injury.Cold hardiness testing was performed over two consecutive weeks (October 18-28, 2021) to accommodate the large number of samples.Test temperatures were −25°C and −30°C while control samples were placed in a fridge at 4°C.A small amount of silver iodide to nucleate ice crystals and 0.2 mL distilled H 2 O were added to samples before freezing.Samples were exposed to test temperatures for 1 h in a programmable temperature chamber, Tenney model T20C-3.After freezing, 3 mL of distilled water was added and samples were shaken for 1 h.Electrical conductivity was measured on test and control samples after freezing, and again after heat killing at 95°C in a laboratory oven overnight, using Amber Science Inc. Model 2052 Digital Electrical Conductivity meters.Cold injury was estimated as the percentage of cellular damage using the ratio of electrolytic leakage after freezing relative to total leakage after heat killing at 95°C between test and control samples.The cold injury damage incurred by each seedling at both test temperatures was calculated as a percentage relative to unfrozen control samples using an index of cold injury (I) (Flint et al. 1967).We averaged the values of I between the two test temperatures and used this mean value for our analyses.

| Population Phenotypes and Variance Components
We used a linear mixed-effect model using residual maximum likelihood to estimate best linear unbiased estimates (BLUEs) of natural populations and selected full-sib families with ASReml-R version 4.0 (Butler et al. 2018).The following mixedeffects model was used: where Y ij is the phenotype corresponding to individual i from the population j and block k; is the phenotypic global mean across all individuals (fixed intercept); P is the coefficient for the fixed effect of the jth population; B is the random effect of the kth block; and ijk is the random residual error term.To account for spatial effects within the experiment, row and column positions within the experiment were incorporated into the random residual term with a correlation structure.Models including the spatially correlated residual term consistently reduced Akaike's information criterion (AIC) values. (1) We estimated population variance components and standard errors for each trait by fitting mixed-effects models.The model form used was the same as Equation ( 1) but with all independent variables set as random effects.

| Population Differentiation
We quantified differentiation among populations for each phenotypic trait with two metrics: V POP , a measure of the amount of phenotypic variance among populations relative to the total phenotypic variance and Q ST , a standardized measure of the amount of genetic variance among populations relative to the total genetic variance.V POP is calculated as the phenotypic variance among populations divided by the sum of the residual variance and the variance among populations.V POP is often considered as a proxy for Q ST when relationships among individuals within populations are not known.However, it is an underestimation of Q ST because in the denominator, the total phenotypic variance is used (i.e., population and residual component), rather than the sum of the among-population component and twice the additive variance component used to calculate Q ST (Liepe et al. 2016).For comparison, we approximated Q ST as the phenotypic variance among populations divided by the sum of twice the additive variance component estimated from selected full-sib families and the phenotypic variance among populations.Q ST estimates were then qualitatively compared with F ST , a standardized measure of neutral genetic differentiation among populations, estimated as global F ST based on a dataset generated from pool-seq exome capture for western larch populations (n = 40 individuals per population) that included over 1.4 million SNPs after filtering for missing values and minor allele frequency (B.Roskilly and B. Lind, unpublished data).Q ST estimates higher than global F ST were considered indicative of local adaptation for a given trait (Whitlock and Guillaume 2009).

| Phenotypic Clines
We used simple linear regressions in R version 4.2.0 (R Development Core Team 2022) to test associations between BLUEs of the population phenotypes and selected environmental variables.Regression coefficients were considered significant at p < 0.05.

| Comparison of Populations and Selected Full-Sib Families Within Breeding Zones
Best linear unbiased estimates (BLUEs) were estimated for seed source type (natural population or full-sib family) by breeding zone using the following mixed-effects model: where Y ijk is the phenotype corresponding to individual i from seed source type j, breeding zone k, block l; is the phenotypic global mean across all individuals, S is the fixed effect of the jth seed source type, Z is the fixed effect of the kth breeding zone, S*Z the fixed effect of the seed source type by breeding zone interaction, B is the random effect of lth block, and ijkl is the random residual error term including the spatial autocorrelation of individual seedling position.We tested for pairwise differences between BLUEs of seed source type phenotypes within breeding zones using pairwise t-tests with the asremlPlus package for ASReml-R version 4.0.

| Trait Correlations
To test for correlations among traits, Pearson correlation coefficients were estimated with BLUEs of the population and family phenotypes.Pearson correlation coefficients were significant at p < 0.05.

| Additive Genetic Variance and Heritability
We estimated additive genetic variance and narrow-sense heritability in selected full-sib families with the individual (animal) model approach using the following mixed-effects model: where Y ijk is the phenotype corresponding to individual i from block j and individual additive genetic value k; is the phenotypic global mean across all individuals, B is the random effect of jth block, G is the random effect of the kth individual additive genetic value derived from the relationship matrix based on the pedigree of the full-sib family structure, and ijk is the random residual error term.
Narrow-sense heritability (h 2 ) for each trait was estimated as: where 2 a is the additive genetic variance estimated from Equation (3) and 2 p is the total phenotypic variance.

| Population Differentiation and Clines
The proportion of phenotypic variance among populations (V POP ) was generally low but varied among traits, ranging from 0.02 to 0.16 (Table 1; Figure S3).The proportion of additive genetic variance among populations (Q ST ) was low to moderate among traits, ranging from 0.10 to 0.28 (Table 1).Q ST estimates were not strongly correlated with V POP due to differences in additive genetic variation between traits.
Population differentiation for timing of bud flush was low (V POP = 0.09, Q ST = 0.20; Table 1) and was weakly associated with longitude (R 2 = 0.14, p = 0.01; Figure 2).Bud flush timing was not correlated with annual growth or date of final bud set, either among populations or among selected full-sib families (Figure S4).The timing of final bud set also had weak population differentiation (V POP = 0.11, Q ST = 0.11; Table 1).Final bud set was moderately correlated with annual growth among populations but not correlated among selected full-sib families (Figure S4).Final bud set was weakly associated with elevation (R 2 = 0.14, p = 0.01; Table 2) and latitude (R 2 = 0.08, p = 0.05; Table 2).
Population differentiation for fall cold injury was very weak (V POP = 0.02, Q ST = 0.10; Table 1).Fall cold injury was not significantly associated with any climatic or geographic variable (Table S2).The range of values for population phenotypes was small, with only 16% average cold injury separating the least and most injured populations (Table 1, Figure S3).

| Comparisons Between Natural Populations and Selected Full-Sib Families
Selected full-sib families from the breeding program averaged 21% greater annual growth than natural populations in both breeding zones (Figure 3).Annual growth had moderately low heritability (h 2 = 0.17; Table 3).Selected full-sib families had greater early and late growth compared to natural populations (Figure S5).Late growth had moderate heritability (h 2 = 0.26) and early growth had low heritability (h 2 = 0.08; Table 3).Bud flush timing had low heritability (h 2 = 0.13; Table 3) and did not differ between selected full-sib families and natural populations (Figure 3).Final bud set was delayed by 4.3 days on average in selected full-sib families compared to natural populations in both breeding zones (Figure 3; Table S3) and had moderately high heritability (h 2 = 0.37; Table 3).Fall cold hardiness did not differ significantly between selected full-sib families and natural populations (Figure 3) and had low heritability (h 2 = 0.07; Table 3).

| Weak Signals of Local Adaptation Among Natural Populations
We studied genetic variation in traits relevant to climate in natural populations and selectively bred families of western larch to better understand population differentiation, the quantitative genetics of climate-relevant traits, and the implications for assisted migration and breeding for future climates.We found low differentiation among populations for all traits, with the proportion of phenotypic variance (V POP ) ranging from 0.02 to 0.16.Estimates of the proportion of additive genetic variance due to population differences (Q ST ) were not strongly correlated with V POP due to differences in additive genetic variation among traits and ranged from 0.10 to 0.28 (Table 1).Although low, Q ST estimates were higher than differentiation among populations for selectively neutral loci that had an estimated global F ST = 0.01 (B.Roskilly and B. Lind, unpublished data), suggesting these traits are weakly locally adapted.Weak clines in height growth, bud flush, and bud set were found with longitude and elevation (Figure 2) but surprisingly not with the climatic variables for source populations that we expected to be drivers of local adaptation in temperate conifers (Table S2).We found that cold hardiness had the weakest population differentiation and had no significant clines.Comparisons between full-sib families and natural populations indicate that selective breeding for faster growth has delayed bud set but maintained trait variation in bud flush phenology and cold hardiness.Collectively, our results indicate weak local adaptation in western larch and suggest phenotypic plasticity may be important in accommodating environmental variation across the species range (Leites and Benito Garzón 2023).
Population differentiation and clines in height growth typically reflect adaptation to growing season conditions (Alberto et al. 2013;Aitken and Bemmels 2016).In most temperate evergreen conifers, height growth tends to be faster in populations from warmer climates, reflected in clines with latitude, elevation or climate (Alberto et al. 2013).In this study, population differentiation for annual growth in western larch populations was low (Table 1), weakly correlated with longitude (Figure 2) and uncorrelated with climatic variables (Table S2).
Comparisons of patterns in this study with the well-studied species Douglas-fir (Pseudotsuga menziesii) are informative as Larix and Pseudotsuga are sister genera in the Pinaceae with sympatric populations in their natural ranges.Population differentiation in height growth in western larch was lower in this study than for populations of Douglas-fir (V POP = 0.17; St Clair  Liepe et al. 2016).Weak population differentiation in height growth for western larch was also found in a provenance trial in northern Idaho (Rehfeldt 1995) and in a series of provenance trials in southern British Columbia (Leites et al. 2012).Clines in height growth with elevation were detected in a northern Idaho provenance trial at 4 years (Rehfeldt 1995) and in 2-year nursery bed trials (Rehfeldt 1982).Though we sampled a broad portion of the species range, it is possible that the less intensive distribution of sampled populations over topographically complex and climatically variable regions of the range or the warmer climate at the study site resulted in weaker associations with elevation than previous studies.These findings serve as an important reminder that signals of local adaptation found in common garden studies can be affected by the climatic variability represented in sampled populations and the test site climate.
As an early-successional, shade-intolerant species, indeterminate free growth is an important trait that enables western larch trees to take advantage of favorable current growing season conditions and gain a dominant or co-dominant position in the canopy.Unlike most conifers in the Pinaceae, free growth extends beyond the juvenile stage and remains a significant component of growth potential in mature western larch trees (Schmidt and Shearer 1990).We could not separate determinate and free growth in this study because preformed and neoformed leaves were not readily distinguishable based on morphological differences, and so instead analyzed early and late growth.Early growth therefore includes both determinate growth from shoot primordia that overwintered in buds and indeterminate free growth from neoformed leaves and internodes during the growing season.Like annual growth, we found early growth had low population differentiation (Table 1) and weak clines with longitude and the length of the frost-free period (Table 2).Faster early growth in eastern populations may reflect an adaptation to compensate for shorter growing seasons.Faster height growth rates in populations with shorter growing seasons have been found in other Pinaceae species, such as Picea sitchensis (Mimura and Aitken 2007).
Late growth included two forms of indeterminate growth: free growth and lammas growth.The latter results from reflushing after bud set in the same growing season.Late growth had weak clines with elevation and latitude following clines in final bud set (Table 2).Populations from lower elevations may have a greater capacity for growth and later bud set to take advantage of favorable late-season growing conditions such as those found at the study site.Earlier bud set in high elevation populations likely reflects adaptation to avoid early fall frost risks and shorter growing seasons.Late growth and bud set also had a weak negative relationship with latitude (Table 2), such that populations from higher latitudes had later bud set.While this was an unexpected pattern, given the limited variation in photoperiod over the latitudinal gradient sampled, it may be due to a moderate correlation between elevation and latitude of provenances (Figure S6), with higher latitude populations growing at lower elevations with less risk of early fall frost (Sang et al. 2021).It is also possible that the favorable late-season growing conditions at the study site reduced population differentiation and weakened the strength of clines in final bud set by promoting reflushing after initial bud set.Expression of lammas growth at the study site was likely much greater than the frequency of lammas growth in their source environments (B.Jaquish, personal communication).However, our results are consistent with previous provenance trials that have generally found elevational clines of only modest slope in western larch, with populations on average needing to be separated by at least 400 m elevation to exhibit significant genetic differentiation (Rehfeldt and Jaquish 2010).
Deciduous tree species tend to break bud earlier than evergreen species because they need regrow their entire canopy to start photosynthesizing at the beginning of each growing season, which puts them at greater risk of late spring frosts (Panchen et al. 2014).We found that bud flush in western larch seedlings was on average 23 days earlier than interior variety Douglasfir seedlings from sympatric populations at the same common garden site in the same year (Nuhu 2022).Bud flush in our western larch seedlings occurred by March 21st, on average, nearly a month earlier than average bud flush in a previous provenance trial within the species range in northern Idaho (Rehfeldt 1992), suggesting that bud flush was substantially advanced by warmer temperatures at the study site.Bud flush was later in eastern populations following a weak cline with longitude (Figure 2).This result is not consistent with the general pattern found in evergreen conifers, where populations from colder or more continental climates typically have earlier bud flush due to lower chilling or heat sum requirements (Howe et al. 2003).Populations from the east may have higher chilling requirements that were met later at the warmer study site, tempering the advancement of bud flush.Though we found only a weak signal of local adaptation in bud flush, managers should be cautious of increasing risks of late spring frosts as a result of assisted migration and climate change given western larch's early bud flush phenology.
The weak population differentiation and lack of clines for cold hardiness in western larch are striking given that cold hardiness typically exhibits the strongest population differentiation and climatic clines of all seedling traits studied in temperate conifer species (Alberto et al. 2013;Leites and Benito Garzón 2023).Even though test temperatures resulted in optimal, intermediate levels of damage expected to maximize phenotypic variation using well established methods, only 16% difference in cold injury separated the least and most cold-injured populations.This is consistent with the lack of significant population differentiation or clines for visible winter cold injury found in situ in a previous provenance trial at a substantially colder test site in northern Idaho (Rehfeldt 1995).Our results suggest western larch may be globally rather than locally adapted to cold temperatures.Even in the warm study site in mid-October, we needed to test branch tissue samples at temperatures of −25 and −30°C to achieve targeted cold injury values of approximately 50%.These temperatures are well below the lowest average minimum temperatures for October at the five coldest source locations (−3.5°C) suggesting a wide safety margin in fall cold hardiness.Global adaptation to cold temperatures in western larch is consistent with the boreal evolutionary history of the Larix genus (LePage and Basinger 1995).

| Selective Breeding for Faster Growth Delayed Bud Set With No Strong Effects on Bud Flush or Cold Hardiness
Selected full-sib families had 21% greater annual height growth on average compared to natural populations (Figure 3), corroborating reported genetic gains in growth in a warmer test site located outside of established breeding zones (Forest Genetics Council of British Columbia 2021).Selection for greater annual growth has  Note: σ 2 a = additive genetic variance estimated from pedigree, σ 2 p = total phenotypic variance for the trait, narrow-sense heritability (h 2 ), and standard error (se) for each phenotypic trait estimated from 28 full-sib families from the breeding program of British Columbia.increased early and late growth in both breeding zones (Figure S5).growth had higher heritability than annual growth and early growth (Table 3).Early growth heritability may have been underestimated in part due to the inability to distinguish predetermined and free growth in study.On the other hand, estimated heritability of late growth in this study was comparable to a previous estimate based on 52 families from three natural populations in the Rocky Mountains (h 2 = 0.20; Rehfeldt 1992).Our heritability estimates for late growth and bud set timing are comparable to those for lammas growth and bud set timing, respectively, in Douglas-fir (Kaya, Adams, and Campbell 1994;Howe et al. 2003).Late growth is increased through a delay of 4.3 days on average in final bud set in selected families compared to natural populations from both breeding zones (Figure 3; Table S3).
Selective breeding has had little impact on bud flush and fall cold hardiness in western larch (Figure 3).Selected families in the eastern zone had slightly elevated cold injury compared to natural populations but the difference was only marginally significant (p = 0.06) and the effect was very small (Figure 3).Cold hardiness and bud flush also had relatively low heritability (Table 3), lower than average cold hardiness heritability estimates for Douglasfir and lower than average estimates for bud flush heritability in other temperate conifer and broadleaf species (Howe et al. 2003).Consistent with these results, bud flush was not correlated with annual growth or final bud set in natural populations or among selected families (Figure S4), likely because indeterminate growth decouples bud flush from annual growth and bud set in western larch.We tested cold hardiness once over a 2-week period in late October during a critical period of cold acclimation in the fall.It is possible that differences between seed sources may be greater at other time periods, for example, earlier in the fall or later in spring.Nonetheless, our findings suggest that selective breeding for faster growth has mostly maintained cold hardiness compared to natural populations, consistent with findings for lodgepole pine and interior spruce (hybrid populations of Picea glauca, Picea engelmannii, and Picea sitchensis; MacLachlan et al. 2017; MacLachlan, Yeaman, and Aitken 2018).However, genomic studies in lodgepole pine indicate that future responses to selection for growth may be constrained by correlated responses among growth, phenology and cold hardiness due to extensive pleiotropy (Maclachlan et al. 2021).

| Potential Causes and Implications of Weak Local Adaptation in Western Larch
The weak signals of local adaptation to climate in western larch we found are consistent with previous work (Rehfeldt 1995;Rehfeldt and Jaquish 2010), though the clines found in our study may differ to the warmer and milder climate of the study site or the distribution of sampled populations.High phenotypic variance within populations for traits in this study (Figure S3) as well as previous work suggest that local adaptation is not limited by low standing genetic variation (Rehfeldt 1992(Rehfeldt , 1995)).Weak signals of local adaptation in this species are likely due to the relatively narrow climatic range inhabited by the species, which could be limited by ecological factors such as disturbance and succession (Rehfeldt and Jaquish 2010).Weak population structure due to population history and few geographic barriers to gene flow across the natural range may also contribute to weak signals of local adaptation.Neutral genetic differentiation among populations was notably low; F ST = 0.01 based on over 1.4 million SNPs (B. Roskilly and B. Lind, unpublished data).This estimate is lower than that estimated with a similar set of SNPs for interior Douglas-fir populations, which are sympatric with western larch populations but span a wider geographic and climatic range (F ST = 0.03; Candido-Ribeiro and Aitken 2024).Western larch also tends to be a minor component of the mixedconifer stands it inhabits, with smaller local population sizes (Schmidt and Shearer 1990).Low density stands may receive relatively high proportions of pollen originating from other locations which would contribute to low differentiation among populations.
This common garden experiment was conducted at a site outside the current natural species range to simulate the effects of warming at a sensitive early developmental stage.Short-term common garden experiments can provide valuable insight into intraspecific genetic variation in climate-relevant traits, facilitating in-depth phenotyping not feasible in remote field locations, and can indicate potential drivers of clinal variation.However, to robustly test how far populations and species can be transferred to track climate change with assisted migration, populations and species should be tested in a wide range of climates in field reciprocal transplant studies (Wadgymar et al. 2022;Leites and Benito Garzón 2023).These studies should be monitored long term, as evidence of maladaptation can take decades to become apparent in long-lived tree species (St. Clair et al. 2020) and could lead to more conservative climate-based seed transfer guidelines (O'Neill et al. 2014).As drought becomes an increasingly dominant driver of tree mortality globally (Hartmann et al. 2022), intraspecific genetic variation and phenotypic plasticity of drought relevant traits also need to be better understood to inform assisted migration and breeding strategies (Moran et al. 2017;Candido-Ribeiro and Aitken 2024).
The weak signals of local adaptation we found among western larch populations suggest limited benefits and risks of assisted gene flow among populations within the species range for mitigating maladaptation with climate change compared with many other tree species.Assisted migration of this species beyond the natural range is already underway based on previous research (Rehfeldt and Jaquish 2010;O'Neill et al. 2017).Our results also suggest that the selection of seed sources for planting in locations outside of the range will be less critical for western larch than for more locally adapted species.We hypothesize that western larch is a species that relies more on phenotypic plasticity than local adaptation to accommodate environmental variation.The ability of this species to respond to climate change may therefore be bounded by the climatic limits of phenotypic plasticity and its interplay with adaptation under changing climates.

FIGURE 1 |
FIGURE 1 | Distribution of 52 sampled natural populations (green circles) and the natural range of western larch in green and location of common garden test site (gold star) in western North America.

FIGURE 2 |
FIGURE 2 | Regressions of bud break and longitude, annual height growth and longitude, late growth and elevation, bud set and elevation among populations (n = 52).Bud flush and final bud set were analyzed as day of year (doy).Shaded areas represent 95% confidence intervals and regressions are significant at a level of p < 0.05.

FIGURE 3 |
FIGURE 3 | Distribution of individual data and estimated population and family phenotypes (bold points) for annual height growth, bud set and cold injury data compared between selected full-sib families (n = 28) and natural populations (n = 23) from the eastern and western breeding zones of British Columbia.Significance of t-tests: *p < 0.05, **p < 0.01, ***p < 0.001.

TABLE 1 |
Population differentiation and variance components.Proportion of phenotypic variance among populations (V POP ), additive genetic variance among populations (Q ST ), and range of estimated population phenotypes (BLUEs).Bud flush and final bud set were analyzed as day of year (doy).

TABLE 2 |
Coefficient of determination (R 2 ) and slope values for significant regressions between phenotypic traits and geographic or climatic variable of source populations at significance level of p < 0.05.

TABLE 3 |
Trait heritability and variance components for selected full-sib families.